

# Code to generate graphs of outcomes under different compliance regimes
# Macartan Humphreys, Raul Sanchez de la Sierra, and Peter van der Windt 

# Parameters used for figure
b		 = .42 # Share fishers
g 		 = .61 # Share new discoveries
q 		 = .5 	# Probability hypothesis is true
p_orig_T = .8
p_orig_F = .05
p_fish_T = .5  # TO be interpreted as the probability of getting a + result by fishing GIVEN that the original approach failed
p_fish_F = .4  # TO be interpreted as the probability of getting a + result by fishing GIVEN that the original approach failed
p_modi_T = .9 # TO be interpreted as the probability of getting a + result by design modification GIVEN that the method was changed
p_modi_F = .05 # TO be interpreted as the probability of getting a + result by design modification  GIVEN that the method was changed 


number.complying = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {(1-b)*(1-g) + b*(q*(p_orig_T +(1-p_orig_T)*(1-p_fish_T))+(1-q)*(p_orig_F +(1-p_orig_F)*(1-p_fish_F)))}


number.pos.complying = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
		(1-b)*(1-g)*(q*p_orig_T + (1-q)*p_orig_F) + b*(q*(p_orig_T)+(1-q)*(p_orig_F ))}

# (number) (H is really T or F) and (result is POS or NEG) and analysis is (COMPLIANT or NOT) 
		
n.T.pos.comp = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	q*p_orig_T*(b + (1-b)*(1-g))	
	}			
n.T.neg.comp = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	q*(1-p_orig_T)*(b*(1-p_fish_T) + (1-b)*(1-g))	
	}
n.F.pos.comp = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	(1-q)*p_orig_F*(b + (1-b)*(1-g))	
	}		
n.F.neg.comp = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	(1-q)*(1-p_orig_F)*(b*(1-p_fish_F) + (1-b)*(1-g))	
	}
number.complying(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)
n.T.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.T.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.F.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)
	
n.T.pos.nonc = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	q*((1-b)*(g*(p_modi_T)) + b*(1-p_orig_T)*p_fish_T)
	}		
n.T.neg.nonc = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	q*((1-b)*(g*(1-p_modi_T)) + b*0)    # Note cheats do not produce noncompliant negative results
	}
n.F.pos.nonc = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	(1-q)*((1-b)*(g*(p_modi_F)) + b*(1-p_orig_F)*p_fish_F)
	}		
n.F.neg.nonc = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	(1-q)*((1-b)*(g*(1-p_modi_F)) + b*0)    # Note cheats do not produce noncompliant negative results
	}
		
# Check	
n.T.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.T.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.F.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.T.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.T.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.F.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
n.F.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)
	
share.compliant.pos= function(b,g) {      # invariant to b,g
	(n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +  n.T.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) )/
	(n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.T.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.F.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.T.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.compliant.pos.false= function(b,g) {      # invariant to b,g
	n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)/
	(n.T.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.compliant.neg.false= function(b,g) {      # the more cheaters thre are  the less likely that the compliant negatives will  be believable; condition for this?   
												# if p_fish_T > p_fish_F, which seems reasonable 	
	n.T.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)/
	(n.T.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.noncompliant.pos= function(b,g) {      # invariant to b,g
	(n.F.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) + n.T.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) )/
	(n.F.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.T.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.F.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.T.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.noncompliant.pos.false= function(b,g) {      # invariant to b,g
	n.F.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)/
	(n.T.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.noncompliant.neg.false= function(b,g) {      # invariant to b,g
	n.T.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)/
	(n.T.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

	
# No registration req: Note assumption is that behavior of types is no different but they do not separate in the same way	
# Also no change in publishers decisions
###############################################################################################

n.T.pos = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	n.T.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.T.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) 	
	}			
n.T.neg = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	n.T.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.T.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) 	
	}
n.F.pos = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	n.F.pos.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.pos.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) 	
	}		
n.F.neg = function(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) {
	n.F.neg.nonc(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.neg.comp(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) 	
	}

	
share.pos= function(b,g) {      # invariant to b,g
	(n.F.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) + n.T.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) )/
	(n.F.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.T.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.F.neg(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	 n.T.neg(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.pos.false= function(b,g) {      # invariant to b,g
	n.F.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)/
	(n.T.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.pos(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}

share.neg.false= function(b,g) {      # the more cheaters thre are  the less likely that the compliant negatives will  be believable; condition for this?   
												# if p_fish_T > p_fish_F, which seems reasonable 	
	n.T.neg(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F)/
	(n.T.neg(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F) +
	n.F.neg(b,g,q, p_orig_T, p_orig_F, p_fish_T, p_fish_F, p_modi_T, p_modi_F))
	}


#  Strict registration
##################################################################################

share.pos.S= function(b,g) {      		# invariant to b,g
	q*p_orig_T +(1-q)*p_orig_F +b*0
	}

share.pos.false.S= function(b,g) {      # invariant to b,g
	(1-q)*p_orig_F/ (q*p_orig_T +(1-q)*p_orig_F) +b*0
	}

share.neg.false.S= function(b,g) {      # the more cheaters thre are  the less likely that the compliant negatives will  be believable; condition for this?   
										# if p_fish_T > p_fish_F, which seems reasonable 	
	q*(1-p_orig_T)/ (q*(1-p_orig_T) +(1-q)*(1-p_orig_F)) +b*0
	}


### Graph Results	
##################################################################################

gl = .1
gh = .9
nlines = 10
grys = 	rev(gray.colors(nlines, end=.8))	
int = (gh-gl)/(nlines-1)

quick = function(f, gl, gh, main=""){
	plot(b, f(b,gh), type="l", main=main, ylim=c(0,1), ann=F, col=grys[1]); title(main=main, font.main=1)
	for(i in 1:(nlines-1)){lines(b, f(b,gh - i*int), col=grys[i+1])}		
	}		
	
png("compliance.png", width=12, height=10,  res = 300, units = "in")
	par(oma=c(2,2,4,2))
	par(mfrow=c(3,4))	
	b = seq(0,1,.05)
	quick(share.pos, gl, gh, main="Probability of positive result")
	quick(share.pos.S, gl, gh, main="Probability of positive result")
	quick(share.compliant.pos, gl, gh, main="Probability of positive result \n among compliers")
	quick(share.noncompliant.pos, gl, gh, main="Probability of positive result \n among noncompliers")
	quick(share.pos.false, gl, gh, main="Probability that a \n positive result is false")
	quick(share.pos.false.S, gl, gh, main="Probability that a \n positive result is false")
	quick(share.compliant.pos.false, gl, gh, main="Probability that a positive result is \n false among compliers")
	quick(share.noncompliant.pos.false, gl, gh, main="Probability that a positive result is \n false among noncompliers")
	quick(share.neg.false, gl, gh, main="Probability that a negative  \n  result is false")
	quick(share.neg.false.S, gl, gh, main="Probability that a negative  \n  result is false")
	quick(share.compliant.neg.false, gl, gh, main="Probability that a negative result \n is false among compliers")
	quick(share.noncompliant.neg.false, gl, gh, main="Probability that a negative result \n is false among noncompliers")
	mtext(expression(paste(beta, " (Share of fishers in population)")), SOUTH<-1, line=0.4, cex=1.2, outer=TRUE)	
	mtext("No Registration", TOP<-3, 			at=.14, line = 1.1, cex=1.1, outer=TRUE)
	mtext("Binding Registration", TOP<-3, 		at=.39, line = 1.1,cex=1.1, outer=TRUE)
	mtext("Nonbinding Registration", TOP<-3, 	at=.64, line = 1.4, cex=1.1, outer=TRUE)
	mtext("Nonbinding Registration", TOP<-3, 	at=.88, line = 1.4,cex=1.1, outer=TRUE)
	mtext("(Compliers)", TOP<-3, 					at=.64, cex=1.1, line = -.1,outer=TRUE)
	mtext("(Noncompliers)", TOP<-3, 				at=.88, cex=1.1, line = -.1, outer=TRUE)
dev.off()


